home *** CD-ROM | disk | FTP | other *** search
/ IRIX 6.2 Development Libraries / SGI IRIX 6.2 Development Libraries.iso / dist / complib.idb / usr / share / catman / p_man / cat3 / complib / dhgeqz.z / dhgeqz
Text File  |  1996-03-14  |  11KB  |  265 lines
  1.  
  2.  
  3.  
  4. DDDDHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))                                                          DDDDHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      DHGEQZ - implement a single-/double-shift version of the QZ method for
  10.      finding the generalized eigenvalues  w(j)=(ALPHAR(j) +
  11.      i*ALPHAI(j))/BETAR(j) of the equation   det( A - w(i) B ) = 0  In
  12.      addition, the pair A,B may be reduced to generalized Schur form
  13.  
  14. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  15.      SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB,
  16.                         ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,
  17.                         INFO )
  18.  
  19.          CHARACTER      COMPQ, COMPZ, JOB
  20.  
  21.          INTEGER        IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
  22.  
  23.          DOUBLE         PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), B(
  24.                         LDB, * ), BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, *
  25.                         )
  26.  
  27. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  28.      DHGEQZ implements a single-/double-shift version of the QZ method for
  29.      finding the generalized eigenvalues B is upper triangular, and A is block
  30.      upper triangular, where the diagonal blocks are either 1-by-1 or 2-by-2,
  31.      the 2-by-2 blocks having complex generalized eigenvalues (see the
  32.      description of the argument JOB.)
  33.  
  34.      If JOB='S', then the pair (A,B) is simultaneously reduced to Schur form
  35.      by applying one orthogonal tranformation (usually called Q) on the left
  36.      and another (usually called Z) on the right.  The 2-by-2 upper-triangular
  37.      diagonal blocks of B corresponding to 2-by-2 blocks of A will be reduced
  38.      to positive diagonal matrices.  (I.e., if A(j+1,j) is non-zero, then
  39.      B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be positive.)
  40.  
  41.      If JOB='E', then at each iteration, the same transformations are
  42.      computed, but they are only applied to those parts of A and B which are
  43.      needed to compute ALPHAR, ALPHAI, and BETAR.
  44.  
  45.      If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal
  46.      transformations used to reduce (A,B) are accumulated into the arrays Q
  47.      and Z s.t.:
  48.  
  49.           Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)*
  50.           Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)*
  51.  
  52.      Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
  53.           Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
  54.           pp. 241--256.
  55.  
  56.  
  57.  
  58.  
  59.  
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. DDDDHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))                                                          DDDDHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))
  71.  
  72.  
  73.  
  74. ARGUMENTS
  75.      JOB     (input) CHARACTER*1
  76.              = 'E': compute only ALPHAR, ALPHAI, and BETA.  A and B will not
  77.              necessarily be put into generalized Schur form.  = 'S': put A and
  78.              B into generalized Schur form, as well as computing ALPHAR,
  79.              ALPHAI, and BETA.
  80.  
  81.      COMPQ   (input) CHARACTER*1
  82.              = 'N': do not modify Q.
  83.              = 'V': multiply the array Q on the right by the transpose of the
  84.              orthogonal tranformation that is applied to the left side of A
  85.              and B to reduce them to Schur form.  = 'I': like COMPQ='V',
  86.              except that Q will be initialized to the identity first.
  87.  
  88.      COMPZ   (input) CHARACTER*1
  89.              = 'N': do not modify Z.
  90.              = 'V': multiply the array Z on the right by the orthogonal
  91.              tranformation that is applied to the right side of A and B to
  92.              reduce them to Schur form.  = 'I': like COMPZ='V', except that Z
  93.              will be initialized to the identity first.
  94.  
  95.      N       (input) INTEGER
  96.              The order of the matrices A, B, Q, and Z.  N >= 0.
  97.  
  98.      ILO     (input) INTEGER
  99.              IHI     (input) INTEGER It is assumed that A is already upper
  100.              triangular in rows and columns 1:ILO-1 and IHI+1:N.  1 <= ILO <=
  101.              IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
  102.  
  103.      A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
  104.              On entry, the N-by-N upper Hessenberg matrix A.  Elements below
  105.              the subdiagonal must be zero.  If JOB='S', then on exit A and B
  106.              will have been simultaneously reduced to generalized Schur form.
  107.              If JOB='E', then on exit A will have been destroyed.  The
  108.              diagonal blocks will be correct, but the off-diagonal portion
  109.              will be meaningless.
  110.  
  111.      LDA     (input) INTEGER
  112.              The leading dimension of the array A.  LDA >= max( 1, N ).
  113.  
  114.      B       (input/output) DOUBLE PRECISION array, dimension (LDB, N)
  115.              On entry, the N-by-N upper triangular matrix B.  Elements below
  116.              the diagonal must be zero.  2-by-2 blocks in B corresponding to
  117.              2-by-2 blocks in A will be reduced to positive diagonal form.
  118.              (I.e., if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and
  119.              B(j,j) and B(j+1,j+1) will be positive.)  If JOB='S', then on
  120.              exit A and B will have been simultaneously reduced to Schur form.
  121.              If JOB='E', then on exit B will have been destroyed.  Elements
  122.              corresponding to diagonal blocks of A will be correct, but the
  123.              off-diagonal portion will be meaningless.
  124.  
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. DDDDHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))                                                          DDDDHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))
  137.  
  138.  
  139.  
  140.      LDB     (input) INTEGER
  141.              The leading dimension of the array B.  LDB >= max( 1, N ).
  142.  
  143.      ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
  144.              ALPHAR(1:N) will be set to real parts of the diagonal elements of
  145.              A that would result from reducing A and B to Schur form and then
  146.              further reducing them both to triangular form using unitary
  147.              transformations s.t. the diagonal of B was non-negative real.
  148.              Thus, if A(j,j) is in a 1-by-1 block (i.e., A(j+1,j)=A(j,j+1)=0),
  149.              then ALPHAR(j)=A(j,j).  Note that the (real or complex) values
  150.              (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the generalized
  151.              eigenvalues of the matrix pencil A - wB.
  152.  
  153.      ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
  154.              ALPHAI(1:N) will be set to imaginary parts of the diagonal
  155.              elements of A that would result from reducing A and B to Schur
  156.              form and then further reducing them both to triangular form using
  157.              unitary transformations s.t. the diagonal of B was non-negative
  158.              real.  Thus, if A(j,j) is in a 1-by-1 block (i.e.,
  159.              A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=0.  Note that the (real or
  160.              complex) values (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are
  161.              the generalized eigenvalues of the matrix pencil A - wB.
  162.  
  163.      BETA    (output) DOUBLE PRECISION array, dimension (N)
  164.              BETA(1:N) will be set to the (real) diagonal elements of B that
  165.              would result from reducing A and B to Schur form and then further
  166.              reducing them both to triangular form using unitary
  167.              transformations s.t. the diagonal of B was non-negative real.
  168.              Thus, if A(j,j) is in a 1-by-1 block (i.e., A(j+1,j)=A(j,j+1)=0),
  169.              then BETA(j)=B(j,j).  Note that the (real or complex) values
  170.              (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the generalized
  171.              eigenvalues of the matrix pencil A - wB.  (Note that BETA(1:N)
  172.              will always be non-negative, and no BETAI is necessary.)
  173.  
  174.      Q       (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
  175.              If COMPQ='N', then Q will not be referenced.  If COMPQ='V' or
  176.              'I', then the transpose of the orthogonal transformations which
  177.              are applied to A and B on the left will be applied to the array Q
  178.              on the right.
  179.  
  180.      LDQ     (input) INTEGER
  181.              The leading dimension of the array Q.  LDQ >= 1.  If COMPQ='V' or
  182.              'I', then LDQ >= N.
  183.  
  184.      Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
  185.              If COMPZ='N', then Z will not be referenced.  If COMPZ='V' or
  186.              'I', then the orthogonal transformations which are applied to A
  187.              and B on the right will be applied to the array Z on the right.
  188.  
  189.      LDZ     (input) INTEGER
  190.              The leading dimension of the array Z.  LDZ >= 1.  If COMPZ='V' or
  191.              'I', then LDZ >= N.
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. DDDDHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))                                                          DDDDHHHHGGGGEEEEQQQQZZZZ((((3333FFFF))))
  203.  
  204.  
  205.  
  206.      WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
  207.              On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
  208.  
  209.      LWORK   (input) INTEGER
  210.              The dimension of the array WORK.  LWORK >= max(1,N).
  211.  
  212.      INFO    (output) INTEGER
  213.              = 0: successful exit
  214.              < 0: if INFO = -i, the i-th argument had an illegal value
  215.              = 1,...,N: the QZ iteration did not converge.  (A,B) is not in
  216.              Schur form, but ALPHAR(i), ALPHAI(i), and BETA(i), i=INFO+1,...,N
  217.              should be correct.  = N+1,...,2*N: the shift calculation failed.
  218.              (A,B) is not in Schur form, but ALPHAR(i), ALPHAI(i), and
  219.              BETA(i), i=INFO-N+1,...,N should be correct.  > 2*N:     various
  220.              "impossible" errors.
  221.  
  222. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  223.      Iteration counters:
  224.  
  225.      JITER  -- counts iterations.
  226.      IITER  -- counts iterations run since ILAST was last
  227.                changed.  This is therefore reset only when a 1-by-1 or
  228.                2-by-2 block deflates off the bottom.
  229.  
  230.  
  231.  
  232.  
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.